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Abstract 

The  main  theme  of  this  research  is  the  application  of  high  order  accurate 
schemes  to  complicated  flow  problems.  The  advantage  of  using  high  order 
schemes  for  long  time  simulations  is  widely  recognized  by  now.  For  prob¬ 
lems  where  fine  details  of  the  flow  field  have  to  be  captured  accurately  the 
use  of  high  accuracy  schemes  is  mandatory.  These  two  classes  of  problems 
encompass  many  of  the  current  problems  in  scientific  computing. 

High  order  accuracy  methods,  finite  difference  ,  finite  elements  or  spectral 
methods  pose  many  algorithmic  problems.  They  have  less  numerical  dissi¬ 
pation  and  therefore  they  are  less  robust.  In  particular  when  the  flow  to  be 
simulated  includes  a  shock  wave,  special  treatment  has  to  be  given.  The  fa¬ 
mous  Gibbs  phenomenon  seems  to  imply  that  the  presence  of  the  shock  wave 
prohibits  the  possibility  of  applying  spectral  methods  for  problems  with  dis¬ 
continuous  solutions.  We  have  shown  that  this  is  the  wrong  interpretation  of 
the  Gibbs  phenomenon  and  applied  the  resolution  of  this  phenomenon  in  this 
research  effort.  Thus  in  this  research  we  present  the  development  and  appli¬ 
cation  of  spectral  shock  capturing  techniques  as  well  as  ENO  finite  difference 
and  finite  element  methods  for  realistic  problems. 

For  genuinely  time  dependent  problems  special  care  has  to  be  taken  in 
order  to  preserve  the  accuracy.  Small  errors  that  do  not  show  up  in  steady 
state  calculations  can  be  amplified  and  ruin  the  total  accuracy.  Time  de¬ 
pendent  boundary  conditions  may  pose  problems.  We  have  addressed  those 
issues  in  our  current  efforts. 

The  recent  advent  of  parallel  computers  poses  a  special  challenge  to  the 
users  of  high  order  methods.  Issues  in  parallel  computing  for  the  simulations 
of  incompressibly  and  compressible  flows  had  been  treated. 
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The  algorithmic  developments  took  place  in  parallel  to  the  applications. 
We  concentrated  on  three  main  applications:  Fuel  air  mixing  enhanced  by 
shock  induced  vortices,  Shock  vortex  interactions  and  flow  past  a  blunt  body. 
This  report  will  be  divided  into  the  following  parts 

1.  Spectral  shock  capturing  techniques,  (p.  3-4) 

2.  High  order  ENO  finite  difference  schemes  for  shock  calculations, 

(p.  5-6) 

3.  Parallel  computing  for  high  order  schemes  (p.  7-10) 

(a)  Spectral  Element  methods  for  incompressible  flows. 

(b)  Spectral  methods  for  compressible  flows 

(c)  ENO  schemes 

4.  Fuel  air  mixing  enhanced  by  shock  induced  vortices,  (p.  11-12) 

5.  Shocks  interacting  with  vortices  (p.  13) 

6.  Compressible  wake  flows  (p.  14-15) 

7.  The  efficient  implementation  of  Spectral  and  high  order  schemes. 

(p.  16) 

8.  Multiscale  computing  (p.  17) 

9.  List  of  publications  resulting  (and  acknowledging)  this  grant,  (p 
I-IV) 

Appendix  A  describes  the  shock  capturing  codes  applied  and  their  imple¬ 
mentation,  using  vector  and  parallel  computers. 
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1.  Spectral  shock  capturing  techniques 


The  problem  in  applying  spectral  methods  to  shock  wave  problems  is  their 
sensitivity  to  discontinuities.  The  presence  of  shock  wave  creates  theoretical 
as  well  as  practical  difficulties. 

Two  major  theoretical  breakthroughs  that  occured  in  the  last  period, 
may  clear  the  way  to  the  successful  implementation  of  spectral  methods  in 
simulating  complicated  shock  waves  interactions. 

The  first  development  involved  the  solution  of  the  Gibbs  phenomenon. 
The  Gibbs  phenomenon  is  related  to  the  well  known  fact  that  the  rate  of 
convergence  of  the  sequence  of  partial  sums  of  the  Fourier  series  representa¬ 
tion  of  a  function  deteriorates  rapidly  when  the  function  has  discontinuities. 
In  essence  the  new  result  says  that  the  first  Fourier  coefficients  of  a  piece- 
wise  analytic  function  (a  finite  number  of  bounded  jump  discontinuities  is 
permitted)  carries  much  more  information  than  can  be  revealed  by  forming 
the  N  th  order  Fourier  partial  sum.  In  fact  it  has  been  shown  that  these  N 
coefficients  contain  enough  information  so  that  one  can  constructively  form 
an  approximation  to  the  unknown  piecewise  analytic  function  which  is  expo¬ 
nentially  accurate  in  N.  This  result  had  been  extended  too  any  polynomial 
method,  such  as  the  commonly  used  Chebyshev  method. 

The  second  development  concerns  the  convergence  of  spectral  methods 
for  non-linear  hyperbolic  equations.  It  has  been  shown  that  with  a  suitable 
addition  of  (spectrally  amall)  artificial  dissipation  the  method  converges. 
Moreover  it  has  been  shown  that  the  usual  filtering  technique  can  be  cast  in 
terms  of  the  above  analized  artficial  dissipation. 

While  the  above  results  create  a  breakthrough  that  indicate  that  high 
order  methods  are  useful  for  shock  wave  calculations,  a  lot  of  work  has  still 
to  be  done  in  making  the  process  efficient.  Stabilizing  the  scheme  in  an 
optimal  way,  and  extracting  the  information  in  an  efficient  way  have  to  be 
further  studied. 

In  a  series  of  papers  ([25]  -  [28])  ,  we  have  addressed  the  issue  of  the  Gibbs 
Phenomenon.  In  an  earlier  paper  we  have  shown  that  an  exponentially  con¬ 
vergent  approximation  in  the  maximum  norm  can  be  reconstructed  from  the 
first  N  Fourier  coefficients  of  a  piecewise  analytic  function.  In  [25]  we  dis¬ 
cussed  resolution  properties  of  this  approximation.  In  [26]  we  have  extended 


3 


the  proof  to  the  case  of  many  discontinuities  and  to  the  case  that  the  Legen¬ 
dre  expansion  coefficients  are  given.  In  [27]  we  discussed  the  case  in  which 
we  know  the  first  N  G eg enbauer  expansion  coefficients  of  a  piecewise  analytic 
function.  In  [28]  we  discusses  the  situation  in  which  we  have  a  good  approx¬ 
imation  to  the  interpolation  polynomial  (or  trigonometrical  polynomial)  of  a 
piecewise  smooth  function. 

In  [39]  It  has  been  shown  that  when  the  spectral  method  is  applied  (with 
the  appropriate  smoothing)  to  a  scalar  conservation  law  the  numerical  solu¬ 
tion  indeed  converges  to  the  entropy  solution.  This  result  was  confirmed  in 
[40]  where  we  have  applied  the  techniques  developed  above  to  extract  expo¬ 
nentially  accurate  information  from  spectral  discretization  of  the  nonlinear 
Burgers  equation. 
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2.  High  order  ENO  schemes 

We  have  been  continuing  our  investigation  of  ENO  schemes  using  point- 
flux  and  TVD  Runge-Kutta  time  discretization  formulations  in  the  following 
directions: 

(1)  Towards  the  convergence  issues,  we  have  investigated  entropy  consis¬ 
tency  for  high  order  Hermite  type  finite  difference  and  discontinuous  Galerkin 
methods  as  a  first  step.  We  have  been  able  to  obtain  [5]  a  cell  entropy  in¬ 
equality  for  the  square  entropy,  for  such  high  order  schemes  (no  restriction 
on  order  of  accuracy)  without  resorting  to  the  help  of  nonlinear  limiters.  The 
result  is  valid  for  multi  space  dimensions  with  arbitrary  triangulations  and 
for  any  fluxes  (no  restriction  on  convexity).  This  is  a  significant  improve¬ 
ment  in  obtaining  cell  entropy  inequalities  since  all  the  previous  work  in  this 
direction  must  either  resort  to  modifications  to  existing  limiters  or  resort  to 
complicated  global  analysis,  and  be  restricted  to  second  order  accuracy.  We 
plan  to  continue  the  investigation  towards  full  convergence  proofs; 

(2)  Application  of  ENO  schemes  to  combustion  problems.  The  first  test 
case  is  shock  interaction  with  hydrogen  bubbles.  Different  configuration  of 
bubbles  is  studied  to  see  the  effect  of  shock  interaction.  High  order  accuracy  is 
crucial  in  this  problem  due  to  the  detailed  structures  of  the  solution  behind 
the  shock.  This  project  is  on  going.  Stiff  source  terms  must  be  treated 
adequately  for  future  tests; 

(3)  Comparison  of  two  different  formulation  of  ENO  schemes:  finite  vol¬ 
ume  vs  finite  difference.  The  comparison  is  performed  on  problems  related 
to  curved  boundary,  inflow  outflow  boundary  conditions,  shocks  oblique  to 
the  grids,  and  CPU  time.  It  is  found  out  that  finite  difference  version  of 
ENO  scheme  (the  one  based  on  point  values  and  numerical  fluxes  of  Shu  and 
Osher)  is  much  cheaper  to  run,  and  can  obtain  results  comparable  to  finite 
volume  ENO  in  most  of  the  test  cases. 

(4)  ENO  schemes  based  on  rational  function  building  blocks  are  being 
investigated.  As  a  first  step  we  are  investigating  approximation  results.  The 
potential  here  is  that  in  most  rapid  transition  regions,  rational  functions  ap¬ 
proximate  the  true  solution  better  than  polynomials.  This  project  is  on  going. 
In  [8]  ,  we  have  performed  an  extensive  comparison  of  the  two  formulations  of 
ENO  schemes:  the  cell-averaged  version  first  developed  by  Harten,  Engquist, 
Osher  and  Chakravarthy,  and  the  point  value  version  first  developed  by  Shu 
and  Osher.  The  results  indicate  that  for  most  test  cases,  the  two  formula- 
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tions  of  ENO  schemes  yields  the  same  accuracy  whereas  the  point  value  ENO 
scheme  is  much  faster  In  [38],  we  have  studied  positivity  preserving  finite  vol¬ 
ume  schemes  in  one  and  two  space  dimensions  for  arbitrary  triangulations. 
The  equations  we  solve  are  Euler  equations  of  compressible  gas,  and  positiv¬ 
ity  is  preserved  for  density  and  pressure.  A  general  framework  and  examples 
are  provided  in  this  paper.  In  [41],  we  have  applied  the  ENO  scheme  to 
the  equations  of  viscoelasticity  with  exponential  fading  memory.  Analytical 
results  about  the  linearized  equation  for  large  time  are  obtained,  and  com¬ 
pared  with  the  numerical  results.  Also,  nonlinear  simulations  for  both  short 
time  and  long  time  are  performed.  The  numerical  method  is  stable  and  pro¬ 
duces  very  good  shock  resolution  for  long  time.  In  [34]  we  have  proved  that 
the  high  order  discontinuous  Galerkin  method,  using  approximate  Riemann 
solvers  satisfies  a  cell  entropy  inequality  for  the  square  entropy,  for  arbitrary 
order  of  accuracy  and  for  arbitrary  triangulations  in  multi  space  dimensions, 
with  or  without  using  the  nonlinear  limiters.  This  compares  sharply  with 
similar  results  for  finite  difference  which  can  be  proven  only  for  much  more 
restricted  cases;  convex,  one  dimension,  and  only  for  special  second  order 
schemes.  Related  stability  issues  and  numerical  aspects  are  also  discussed  in 
[35]. 
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3.  Parallel  computing  for  high  order  schemes 


(a)  Spectral  Elements  for  incompressible  flows 

The  research  objective  here  is  to  develop  parallel  algorithms  for  compu¬ 
tational  fluid  dynamics  (CFD)  which  will  permit  solution  of  incompressible 
flows  with  the  accuracy  and  resolution  demanded  by  large  eddy  simulations 
(LES)  of  turbulent  flows  in  complex  geometries.  The  work  is  subdivided  into 
three  principal  areas:  high-order  flexible  discretizations,  fast  multi-level  it¬ 
erative  solvers,  and  implementation  of  LES  modeling  technology.  All  of  this 
work  is  being  developed  within  the  computational  framework  of  distributed 
memory  architectures  which  provide  a  favorable  price/performance  ratio  for 
this  class  of  problems. 

Our  numerical  approach  is  based  upon  the  spectral  element  method, 
which  retains  many  of  the  essential  features  of  global  spectral  discretizations, 
namely,  rapid  (exponential)  convergence,  minimal  numerical  dissipation  and 
dispersion  per  degree-of-freedom,  and  efficient  tensor  product  factorization 
of  spatial  derivative  operators.  The  computational  domain  is  subdivided 
into  large  macro- elements  and  the  solution,  data,  and  geometry  within  each 
element  are  expressed  in  terms  of  high-order  polynomials.  The  use  of  a 
weighted-residual  procedure  permits  a  reduction  in  inter-element  continuity 
requirements  from  to  (7°,  which  in  turn  leads  to  a  reduction  in  inter¬ 
processor  communication.  The  locally-structured/globally-unstructured,  ap¬ 
proach  of  the  spectral  element  discretization  is  ideally  suited  to  the  two-level 
memory  hierarchy  associated  with  distributed  memory  parallel  computers. 

The  performance  of  general  geometry  incompressible  codes  is  largely  tied 
to  the  speed  of  the  elliptic  solver  for  the  pressure;  development  of  fast  multi¬ 
level  iterative  solvers  is  a  major  focus  of  our  current  efforts.  Presently,  we 
employ  a  two- level  conjugate  gradient  based  solver  which  uses  deflation  to 
project  out  the  coarse  grid  modes,  thereby  reducing  the  condition  number  of 
the  underlying  iteration  matrix.  The  system  is  preconditioned  by  local  finite- 
element  based  operators  which  are  significantly  less  expensive  to  invert  than 
their  high-order  counterparts.  In  effect,  one  can  obtain  high-order  accuracy 
at  low-order  cost.  We  have  addressed  all  of  the  parallel  issues  associated 
with  this  approach,  including  the  use  of  an  efficient  parallel  coarse-grid  solve. 
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and  are  able  to  compute  Navier-Stokes  solutions  at  a  resolution  of  3  million 
gridpoints  on  the  512  node  Intel  Delta  at  the  rate  of  4  minutes  per  time  step 
(6.5  GFLOPS).  To  improve  upon  this  result,  we  are  currently  developing  a 
preconditioner  based  upon  a  Schur- complement  formulation  for  the  interface 
variables.  This  approach  leads  not  only  to  better  conditioning,  but  also 
to  reduced  degrees-of-freedom,  which  in  turn  permits  the  use  of  projection 
techniques  in  which  very  good  initial  guesses  can  be  generated  by  projecting 
the  residual  onto  results  from  previous  time  steps. 

Accurate  numerical  simulation  of  many  high  Reynolds  number  engineer¬ 
ing  flows  will  continue  to  be  limited  by  resolution  requirements  for  the  fore¬ 
seeable  future.  However,  progress  is  being  made  to  the  point  where  the 
combination  of  high- resolution  (  >  10®  gridpoints  )  and  advanced  large-eddy 
simulation  (LES)  models  will  be  able  to  conquer  many  important  problems 
in  the  near  future.  Our  goal  is  to  couple  the  latest  LES  technology  developed 
within  the  turbulence  community  with  a  general  geometry  Navier-Stokes  code 
capable  of  solving  these  demanding  problems. 

Because  of  the  large  number  of  degrees-of-freedom  involved,  effective  use 
of  high-performance  distributed  memory  parallel  architectures  is  essential 
to  economic  resolution  of  these  problems.  Principal  parallel  issues  to  be 
addressed  include:  development  of  coarse-grid  solve  strategies  which  will 
remain  competitive  as  the  number  of  processors  and  dimension  of  the  coarse- 
grid  system  continue  to  increase;  and  development  of  optimal  communication 
strategies  for  the  complex  subdomain  interfaces  arising  from  nonconforming 
discretizations. 


(b)  Spectral  methods  for  compressible  flows. 

Spectral  methods  involve  the  approximation  of  the  unknown  solution  in 
terms  of  global  polynomials.  This  fact  make  them  difficult  to  implement  on 
parallel  computers.  A  popular  method  to  overcome  the  limitations  of  spectral 
methods  is  to  use  multidomain  techniques,  in  which  a  complex  domain  is 
decomposed  into  several  simpler  subdomains.  This  method  has  been  applied 
successfully  to  incompressible  flows  (the  Specral-Element  technique)  or  to 
problems  in  structural  mechanics  (the  h-p  method). 

Multidomain  spectral  methods  are  suitable  for  coarse  grain  parallel  com¬ 
puting,  each  domain  is  assigned  to  a  different  processor. 
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The  main  question  is:  Are  multidomain  methods  efficient?  This  question 
has  not  been  yet  answered  for  those  methods  applied  to  hyperbolic  equations. 
If  we  denote  by  W{p,N)  the  work  involved  in  approximating  k  waves  using 
p  sub  domains  (and  N  points  in  each  domain)  to  obtain  an  error  of  at  most 
e~',  then  W{p,N)  is  minimized  if 


Thus  the  optimal  number  of  subdomain  increases  with  the  complexity  of 
the  problem  (or  number  of  waves)  but  decreases  if  the  required  accuracy 
increases.  This  result  may  serve  as  a  guideline  to  the  optimal  number  of  sub- 
domain.  The  formula  above  can  be  suitably  modified  for  parallel  computers. 

A  key  issue  in  the  application  of  multidomain  spectral  methods  to  the 
numerical  solution  of  hyperbolic  equations  is  the  interface  boundary  condi¬ 
tions.  This  leads  to  the  question  of  the  imposition  of  boundary  conditions, 
both  analytic  and  numerical,  in  the  numerical  solutions  of  systems  of  hyper¬ 
bolic  equations.  For  truly  time  dependent  problems  stability  in  the  classical 
sense  (Lax  and  G-K-S  stability)  is  not  enough.  Even  stable  schemes  may 
exhibit  a  non-physical  growth  in  time.  From  a  practical  point  of  view,  in 
order  to  achieve  reasonable  accuracy  for  large  time,  meshes  much  too  fine  for 
the  computers  available  in  the  foreseeable  future  are  required.  In  the  past 
methods  that  preserve  summation-by-parts  property  of  the  differential  equa¬ 
tion  were  constructed.  Recent  attempts  to  utilize  these  boundary  closures 
to  numerically  solve  a  2  x  2  hyperbolic  system  have  shown  that,  in  certain 
cases,  an  unwarranted  growth  in  time  still  results. 

In  [23],  [7]  we  have  outlined  a  systematic  procedure  for  designing  time- 
stable,  as  well  as  G-K-S  stable  schemes  of  high-order  accuracy.  The  new 
schemes  are  guaranteed  to  be  time-stable  for  any  hyperbolic  system  (as  long 
as  the  system  has  a  bounded  energy).  We  have  extended  this  methodology 
to  Navier  Stokes  equations  in  three  space  dimensions.  We  have  showed  that 
the  SAT  boundary  condition  assures  the  correct  behavior  as  the  Reynolds 
number  tends  to  infinity. 

We  have  carried  out  ([31],  [30])  experiments  with  the  new  boundary  pro¬ 
cedure  indicate  that  there  are  very  suitable  for  parallel  computing.  Indeed 
if  one  assigns  each  domain  to  a  different  processor  then  the  work  can  be 
done  completely  in  parallel.  The  communication  issue  is  related  directly  to 
the  interface  boundary  conditions.  With  the  SAT  procedure  it  seems  that 
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communication  time  may  be  reduced.  This  will  be  one  of  the 

We  also  a  examined  a  method  to  approximate  the  interface  conditions  for 
Chebyshev  polynomial  approximations  to  the  solutions  of  parabolic  prob¬ 
lems,  and  a  smoothing  technique  is  used  to  calculate  the  interface  conditions 
for  a  domain  decomposition  method.  The  method  uses  a  polynomial  of  one 
less  degree  than  the  full  approximation  to  calculate  the  first  derivative  so 
that  interface  values  can  be  calculated  by  using  only  the  adjacent  subdo¬ 
mains.  Theoretical  results  are  given  for  the  consistency  of  the  scheme  and 
practical  results  are  presented.  Computational  results  are  given  for  a  fourth 
order  Runge-Kutta  method  in  two  dimensions  and  for  an  explicit/explicit 
scheme  in  both  one  and  two  dimensions. 


(c)  Parallel  implementation  of  ENO  schemes  on  CM-5 

The  main  cost  in  ENO  schemes  is  in  its  logic  step  in  choosing  local  stencils 
by  comparing  divided  difference  tables  of  the  function.  Although  great  effort 
has  been  spent  on  efficiently  vectorizing  this  part  for  CRAY  supercomputers, 
due  to  the  inevitable  gathering-scattering  process,  ENO  schemes  still  do  not 
run  very  fast  on  CRAY  computers.  Recently  we  have  been  exploring  the 
structure  of  the  ENO  algorithm  to  suit  the  parallel  structure  of  CM-5.  The 
algorithm  has  been  slightly  re-formulated  (in  a  mathematically  equivalent 
way)  to  reduce  communications  and  to  eliminate  communications  between 
other  than  next  neighbors,  at  the  price  of  a  slightly  increased  operation 
count.  Our  CM-5  two  dimensional  ENO  code  for  compressible  Euler  and 
Navier-Stokes  equations  is  a  magnitude  faster  than  our  ENO  code  on  CRAY 
for  a  400x400  grid.  Three  dimensional  code  also  shows  a  speed  up,  although  it 
has  not  been  optimized  yet  (for  three  dimensions,  storage  is  a  big  restriction, 
and  the  structure  of  the  program  must  be  modified  accordingly).  Currently 
we  are  trying  to  improve  the  CM-5  code  and  applying  it  to  reactive  flows. 
This  research  is  reviewed  in  [19],  see  also  Appendix  A. 
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4.  Fuel  air  mixing  enhanced  by  shock  induced  vortices 


In  designing  supersonic  combuster  for  the  next  generation  of  supersonic 
transport,  we  look  for  an  efficient  combuster  such  that  : 

•  allow  better  load  to  weight  ratio  by  carrying  less  fuel, 

•  reduce  chemical  product  that  contribute  to  pollution  due  to  incomplete 
burning  of  fuel. 

One  of  the  technique  currently  under  study  in  enhancing  mixing  of  a  hy¬ 
drogen  jet  (fuel)  and  air  (oxidizer)  is  to  allow  the  existing  shock  inside  the 
combustion  chamber  to  interact  with  the  hydrogen  jet .  By  doing  so,  vorticity 
are  generated  according  to  the  vorticity  equation  of  the  Navier  Stokes  equa¬ 
tion,  The  pressure  gradient  of  the  shock  and  the  density  gradient  between 
the  air  and  hydrogen  provides  an  efficient  mechanism  for  vorticity  genera¬ 
tion  along  the  surface  of  the  hydrogen  jet.  The  vorticity  forced  the  jet  to 
curl  up  and  stressed.  The  increased  surface  area  allows  greater  mixing  of  air 
and  hydrogen  where  combustion  can  take  place.  Even  though  the  problem  is 
three  dimensional  steady  state,  it  had  been  argued  that  the  three  dimensional 
steady  flow  can  be  directly  associated  with  a  corresponding  two  dimensional 
unsteady  calculation  in  this  particular  physical  setting. 

Numerical  simulations  are  a  necessary  part  of  this  investigation.  For  these 
type  of  problem,  it  is  important  to  capture  the  complex  physics  with  high 
accuracy.  Finite  difference  scheme  ,  with  their  inherent  dissipative  nature  for 
stability  reason,  can  only  yields  a  quantitative  result  of  the  flow  fields.  Often, 
for  long  time  integration,  the  loss  of  accuracy  in  the  earlier  stage  affect  the 
development  of  the  flow  in  later  time. 

In  Appendix  A  we  describe  two  codes  for  simulating  the  above  problem. 
A  spectral  shock  capturing  code  and  a  finite  difference  ENO  code.  We  have 
employed  all  the  theoretical  techniques  developed  in  order  to  run  succesfully 
the  codes. 

we  computed  the  solution  of  these  problem  with  various  configuration  of 
multiple  jets  placement.  Our  preliminary  results  found  that  different  ini¬ 
tial  placement  of  the  jets  yields  distinct  final  configurations  (Moreover,  the 
spectral  code  is  fully  capable  of  capturing  the  fine  scale  structure  of  the  in¬ 
teractions,  including  some  that  had  not  been  seen  in  finite  difference  code. 
One  distinct  feature  of  this  calculation  is  the  penetration  of  an  air  stream 
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(heavy  fluid)  into  the  hydrogen  (light  fluid)  causing  instability  that  exhibit 
a  mushroom  shape  structure  (See  Appendix  A).  This  study  can  be  used  to 
guided  researcher  in  develop  fuel  jet  configuration  for  the  scramjet  engine 
with  greater  confident.  Our  goal  in  this  research  is  to  extend  this  problem  to 
a  three  dimensional  steady  and/or  unsteady  simulation  with  full  chemistry 
model. 

A  detailed  report  [19]  describes  the  codes  in  detail. 
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5.  Shock  interacting  with  vortices 


In  designing  supersonic  scramjet  engine  nozzle,  acoustic  radiation  pat¬ 
tern  (sound  wave)  will  be  generated  by  the  interaction  of  shock  and  vortices 
according  to  the  linear  analysis  done  by  Ribner  and  Moore  .  The  sound  wave 
will  have  a  significant  impact  on  the  design  of  engine  that  operated  at  super¬ 
critical  nozzle  pressure  ratios.  Moreover,  sound  wave  (noise)  generated  by 
supersonic  aircraft  should  be  minimized  for  environmental  reason,  as  it  travel 
over  populated  land  mass.  Hence,  a  better  understanding  of  the  mechanism 
of  the  shock  vortices  interaction  is  important. 

The  vortex  is  defined  as  follow  :  the  tangential  velocity  profile  of  the 
vortex  centered  at  (xc,?/c)  in  the  polar  coordinate  is 


C/(r) 


rr(ro  ^  0  <  r  <  Tq  <  ri 

<  rr(r“^  —  tq  <  r  <  ri 

0  r  >  rj 


where  Tq  =  0.2  and  ri  =  1.0  unless  specified  otherwise.  This  vortex  is  rotating 
in  a  counter-clockwise  direction.  Hence,  the  velocity  field  upstream  of  the 
shock  becomes 


u  =  Ur  — U{r)  sin  0 
V  =  Vr  U{r)  cos  6 


where  6  =  tan  ^  f • 

We  used  Chebyshev  collocation  methods  to  solve  the  two  dimensional 
Euler  gas  dynamics  equation  [14].  Using  some  standard  numerical  techniques, 
the  spectral  code  remains  stable  and  highly  efficient,  we  can  computed  the 
solution  in  less  than  10  CPU  minutes  on  a  Cray  2  supercomputer  with  grid 
size  up  to  256x256  Chebyshev  collocation  points.  As  a  test  case,  a  mach  3 
shock  was  propagated  through  air  and  impacted  with  a  vortex  of  strength 
r  =  0.25.  After  reconstruction  of  the  raw  data,  we  are  able  to  recover  the 
solution  that  depicted  strong  acoustic  wave  propagated  downstream  with 
good  accuracy  when  compared  with  ENO  third  order  scheme. 
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6.  Compressible  wake  flows 


We  are  interested  here  in  numerical  simulation  of  a  low  Mach  num¬ 
ber,  compressible,  viscous  flow  past  a  circular  cylinder  at  a  moderately  low 
Reynolds  number.  The  unsteady  wake  generated  by  the  cylinder  has  been 
of  great  interest  to  computational  fluid  dynamicisits  as  well  as  to  theoretical 
and  experimental  aerodynamicits.  The  Reynolds  number  range  between  40 
and  1000  has  been  of  particular  interest  because  it  spans  the  transition  from 
steady  flow  to  an  unsteady  wake  flow  dominated  by  the  period  shedding  of 
vortices  from  the  cylinder.  The  shedding  frequency  of  this  vortices  is  known 
theoretically  and  all  numerical  methods,  applied  to  this  problem,  obtain  this 
frequency  within  reasonable  accuracy.  However,  Srinivassan  measured  (in 
a  wind  tunnel  experiment)  more  than  one  distinct  frequency  in  the  shed¬ 
ding  regime.  In  addition  to  the  vortex  shedding  frequency  he  found  clearly, 
discernible  lower  frequencies  and  concluded  that  this  was  a  feature  of  the 
initial  stage  of  transition  to  turbulence.  A  numerical  simulation,  using  the 
second  order  accurate  MacCormack  scheme  did  find  a  secondary  frequency, 
very  nearly  the  same  found  in  the  wind-tunnel  experiment. 

We  have  developed  at  Brown  the  most  accurate  code  for  this  configura¬ 
tion.  The  code  simulates  the  axial-symmetric  time  dependent  ,  compressible 
viscous  Navier-Stokes  equations.  It  uses  Fourier  methods  in  the  9  direction 
and  Chebyshev  methods  in  the  radial  direction.  This  code  found  no  sec¬ 
ondary  frequencies.  Spectral  methods  are  ,  of  course,  methods  with  high 
order  accuracy,  and  as  any  other  high  order  method,  are  very  sensitive  to  the 
kind  of  outflow  boundary  conditions  treatment.  In  particular  they  tolerate 
only  those  treatments  based  on  the  characteristic  variables.  Low  order  meth¬ 
ods  are  more  tolerant  and  allow  high  variety  of  boundary  conditions.  When 
the  far  field  boundary  conditions  from  the  spectral  code,  were  incorporated 
into  the  second  order  code,  the  secondary  frequency  disappeared.  The  spuri¬ 
ous  frequency  was  generated  by  a  boundary  condition  that  was  tolerated  by 
the  low  order  (and  robust)  MacCormack  scheme  but  not  by  the  high  order 
scheme  used.  It  is  interesting  to  note  that  when  the  wind  tunnel  experiment 
was  repeated  in  a  bigger  wind  tunnel  where  effect  from  the  walls  were  elimi¬ 
nated,  the  secondary  frequency  disappeared.  The  above  example  illustrates 
the  type  of  surprises  one  might  encounter  when  using  low  order  schemes  for 
long  time  integrations.  Stability,  in  the  classical  sense  of  Lax  ,  is  not  enough 
anymore,  one  has  to  be  careful  not  to  have  spurious  phenomena.  It  also  illus- 
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trates  the  fact  that  the  gross  features  of  the  flow  (in  this  case  the  shedding 
frequency)  were  obtained  even  with  low  order  schemes.  It  is  only  when  some 
delicate  features  of  the  flow  were  sought  for  that  an  unphysical  solution  had 
been  observed.  The  fact  that  the  lower  order  scheme  produced  a  spurious 
solution  does  not  contradict  the  fact  that  the  scheme  is  stable  and  converges 
as  the  number  of  grid-points  increases.  The  spurious  solution  will  eventually 
disappear  if  further  grid  refinements  will  be  carried  out.  However  for  real¬ 
istic  grids  it  is  impossible  to  distinguish  between  a  physical  and  a  spurious 
solutions,  a  high  order  scheme  is  less  likely  to  produce  a  spurious  solution. 
With  this  in  mind  we  further  investigated  delicate  physical  questions. 

The  effect  of  heating  of  the  cylinder  on  the  shedding. 

This  is  a  fundamental  question  and  there  are  only  few  wind  tunnel  exper¬ 
iments  (done  by  Srinivassan)  to  indicate  the  possible  effects  on  the  shedding. 

The  interaction  of  acoustic  waves  with  the  flow  and  its  effect  on 
the  shedding. 

A  compact  4-th  order  and  sixth  order  finite  difference  scheme  has  also 
been  incorporated  into  the  spectral  code.  (See  for  a  description  of  the  code 
in  [20]  A  detailed  analysis  of  the  effects  of  the  heating  of  the  cylinder  on 
the  shedding  frequency  had  been  studied  using  the  above  schemes.  It  has 
been  found  that  the  shedding  frequency  decreases  when  the  wire  was  heated. 
Experimental  work  carried  out  at  NASA  Langley  confirmed  the  numerical 
results. 

We  have  also  applied  spectral  multidomain  techniques  to  the  same  geom¬ 
etry  ([30])  in  order  to  be  able  to  port  the  code  easily  to  parallel  computers. 

We  are  now  in  the  process  of  studying  interactions  of  acoustic  waves  and 
heated  cylinders.  The  physical  problem  is  of  great  interest  to  experimentalist 
that  have  to  calibrate  their  wind  tunnel  experiments.  Since  the  probes  are 
heated  cylinders,  the  effect  of  the  acoustic  waves  on  them  has  to  be  taken 
into  account. 
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7.  The  efficient  implementation  of  spectral  and  high  order  schemes 


High  order  finite  difference  methods  were  considered  as  an  alternative  to 
Spectral  methods.  In  [5]  a  series  of  compact  fourth-  and  sixth-order  schemes 
were  developed  and  their  stability  for  mixed  initial-boundary  value  problems 
had  been  verified  using  the  GKS  theory.  In  [6]  The  problem  of  time  stability 
was  also  addressed  and  fourth  order  schemes  that  satisfy  evergy  norm  were 
constructed.  In  [7]  we  had  pointed  out  that  time  stability  for  system  does 
not  follow  from  scalar  time  stability  and  presented  a  systematic  way  for 
constructing  boundary  conditions  for  compact  (and  spectral)  schemes.  This 
is  the  first  work  that  the  issue  of  time  stability  for  systems  has  been  addressed. 

In  [7]  we  have  shown  that  the  conventional  method  of  imposing  time 
dependent  boundary  conditions  for  Runge-Kutte  time  advancement  reduces 
the  formal  accuracy  of  the  space  time  method  to  second  order.  This  counter 
intuitive  result  was  analyzed  and  a  remedy  was  given  for  linear  problems.  A 
partial  remedy  had  been  given  for  some  nonlinear  problems  in  [37]  and  [1].. 

A  step  between  the  theory  and  the  applications  is  the  issue  of  efficient  im¬ 
plementation  of  the  theory.  Progress  had  been  made  in  different  directions. 
In  [18]  the  authors  studied  a  new  method  in  reducing  the  roundoff  error  in 
the  computations  of  derivatives  by  Chebyshev  methods.  In  [17]  mappings 
were  considered  for  Legendre  methods.  In  [15]  we  constructed  a  new  collo¬ 
cation  method  that  uses  Chebyshev  points,  convenient  for  the  applications, 
but  still  has  all  the  theoretical  advanatges  of  the  Legendre  methods.  In 
[31]  we  presented  a  time  stable  open  boundary  conditions  for  the  numerical 
approximation  of  the  Compressible  Navier-Stokes  equations  in  3-D  and  dis¬ 
cussed  their  implementation  within  spectral  and  high  order  finite  difference 
methods. 
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8.  Multiscale  Computing 


Many  relevant  physical  phenomena  involve  infinite  number  of  scales.  In 
nonlinear  problems  it  is  desirable  to  find  a  way  to  take  into  account  the  effect 
of  the  scales  which  are  neglected  on  those  which  are  taken  into  account. 

Two  approaches  have  been  investigated  for  this  type  of  problems:  Non¬ 
linear  Galerkin  Methods  and  Wavelet  based  Schemes. 

In  [29]  we  have  studied  the  implementation  of  the  Nonlinear  Galerkin  in 
the  context  of  collocation  discretizations.  We  have  found  interesting  charac¬ 
terizations,  in  the  physical  space,  of  a  small  scale  function  and  a  large  scale 
function. 

In  [13]  we  presented  an  efficient  pseudospecral  NLG  scheme  for  the  peri¬ 
odic  Burgers  equation,  The  case  of  Chebyshev  approximations  to  nonperiodic 
problems,  in  which  the  concept  of  large  and  small  scales  have  to  be  redefined 
is  currently  under  investigation. 

In  a  thesis  by  L.  Jameson  and  in  a  subsequent  paper  [32],  the  Daubechies 
wavelet  based  differentiation  matrix  was  constructed.  The  relationship  be¬ 
tween  this  matrix  and  finite  difference  methods  was  clarified.  This  serves 
as  a  basis  for  current  work  by  Bruce  Bauer  a  doctoral  student  on  wavelet 
optimized  finite  difference  methods. 

Jameson  [33]  also  analyzed  the  differentiation  matrix  based  on  the  com¬ 
pactly  supported  Daubechies  wavelets.  He  showed  that  in  this  case  there  is 
a  loss  of  the  superconvergence.  The  same  result  holds  (see  [24])  for  the  finite 
element  schemes. 
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Appendix  A 


1  Introduction 


This  research  involves  study  of  combustion  and  mixing  induced  by  the  interaction  of  a  shock 
in  air  with  a  hydrogen  jets.  This  model  problem  is  of  interest  for  the  design  of  air-breathing 
scramjet  engines,  as  it  provides  a  mechanism  for  inducing  millisecond  combustion  times. 
Extremely  rapid  combustion  is  a  necessity  in  an  engine  where  the  reactants  are  escaping 
from  the  reaction  chamber  at  supersonic  speeds. 

In  a  pre- mixed  combustion  process,  the  combustion  rate  is  limited  by  the  rate  at  which 
reactants  diffuse  across  the  air-fuel  interface.  This  diffusion  rate  is  determined  by  the  length 
of  the  interface  and  the  concentration  gradient  across  it.  One  way  to  increase  this  rate  of 
diffusion  is  to  stretch  out  the  interface  by  the  motion  induced  from  a  point  vortex  [2].  At 
the  same  time  this  would  steepen  the  concentration  gradient  across  the  interface  by  bringing 
fresh  reactants  into  contact  with  it.  Marble  [3]  suggests  a  way  to  induce  this  vorticity — to 
have  a  shock  pass  though  a  jet  of  fuel  nearly  perpendicular  to  its  axis.  While  passing  through 
the  region  of  inhomogeneous  density  the  shock  will  produce  vorticity  via  the  term  in 

the  vorticity  equation: 


=  (a; .  V)y  -  cuV  •  1/ + 


Vp  X  Vp 


The  gradient  in  pressure  across  the  shock  in  conjunction  with  gradient  in  fluid  density 
between  the  air  and  hydrogen  produce  a  large  increase  in  vorticity  as  the  shock  passes 
through  the  hydrogen  jet. 

Another  factor  affecting  mixing  in  this  model  system  is  the  fluid  instability  produced  as 
the  denser  air  and  lighter  hydrogen  are  accelerated  by  the  shock  passing  through.  This  is 
a  similar  situation  to  a  heavy  fluid  lying  above  a  lighter  one  in  a  gravitational  acceleration. 
The  heavier  fluid  tends  to  form  fingers  descending  down  into  the  lighter  fluid  below,  due 
to  the  Rayleigh-Taylor  instability.  A  similar  instability,  the  shock-induced  Rayleigh-Taylor 
instability  or  Richtmyer-Meschkov  instability,  results  in  a  similar  fingering  of  air  into  the 
hydrogen  cylinder  as  the  shock  passes  through. 

For  this  type  of  problem,  it  is  important  to  capture  the  complex  physics  with  high 
accuracy.  Accurate  calculations  in  search  of  quantitative  information  about  the  mixing 
require  high  resolution  and/or  high  order  schemes.  Previous  researchers  failed  to  capture 
many  important  features  of  this  flow  due  to  the 


•  inability  of  their  codes  to  adequately  resolve  features  in  the  flow  [4]  [5]  [6]. 

•  inherently  dissipative  nature  of  the  finite  difference  scheme  (FCT  [1,  5],  for  example), 
which  ensured  that  the  computed  results  were  only  qualitatively  accurate.  Often,  for 
long  time  integration,  the  loss  of  accuracy  in  earlier  stages  inhibited  the  development 
of  fine  scale  features  at  later  times. 


•  i  k 


According  to  Marble  et  al.,  the  3-D  steady  shock  and  profile  can  be  simulated  by  a  2-D 
unsteady  shock  provided  several  conditions  can  be  met.  To  make  any  direct  association  of 
the  2-D  unsteady  flow  with  the  3-D  steady  flow  in  such  a  geometry,  it  is  required  that 


•  Pressure  and  Density  jumps  are  matched  across  the  shock. 

•  The  velocity  of  the  feature  {dxldt\2D)  in  2-D  should  be  related  to  the  corresponding 
slope  (dz/dxlso)  in  3-D  linearly 


Idx  dz 

-.tM = "‘s 


ZD") 


where  m  is  a  conversion  factor  and  c  is  the  sound  speed. 


•  To  be  consistent  with  the  pressure  and  density  jump  across  the  shock,  the  motion  of 
the  shock  should  be  matched  cis  well,  hence. 


M2D  =  TOtan(/ii  —  ^). 

where  ^  is  determined  by  the  knowing  Mzd  and  the  turning  angle  8. 


Once  the  conversion  factor  m  is  found  for  a  given  Mzd-,M2Di8^  one  can  compare  the  2-D 
trajectory  of  the  center  of  mass  fraction  from  the  numerical  simulation  with  the  3-D  jet  lift¬ 
off.  They  find  a  good  agreement  for  Mzd  =  6,^  =  tan“^(l/12)  and  M2D  =  1-346,  m  =  9.34. 
This  good  agreement  indicates  that  the  trajectory  of  the  feature  in  2-D  corresponds  well 
with  the  3-D  steady  case  with  the  matching  of  geometry  and  shock  condition. 

We  solved  the  two  dimensional  compressible  Navier  Stokes  equations  with  three  addi¬ 
tional  equations  for  the  conservation  of  mass  of  species  in  order  to  take  into  account  of  the 
production  of  species  by  the  combustion  process.  A  single  step  reversible  chemistry  reaction 
model,  namely,  2H2  -k  O2  <=>  2H2O  is  used.  The  Soret  and  Dufour  effect,  heat  radiant 
effect,  pressure  gradient  diffusion  and  body  force  are  neglected.  The  binary  coefficient  of  all 
species  is  assumed  to  be  equal  and  Lewis  number  Te  =  1.  The  mixture  viscosity  is  deter¬ 
mined  from  the  Wilke’s  law.  The  forward  reaction  rate  of  the  reaction  is  given  by  a  modified 
Arrhenius  law. 

Two  different  numerical  schemes  are  used,  namely,  a  spectral  scheme  [8]  and  the  ENO 
schemes  of  Shu  and  Osher  [7].  We  begin  first  by  discussing  the  spectral  code  in  next  section. 


2  Spectral  code 

•  The  physical  domain  of  the  model  is  x  G  [0  cm,  12.5  cm].y  G  [—6.5  cm,6.b  cm].  The 
symmetry  property  of  the  flow  is  used  for  the  single  jet  configuration.  This  would 
improve  the  efficiency  for  these  cases. 

•  Initially,  a  Mach  2  noi'mal  air  shock  located  at  x  =  5  cm  and  a  hydrogen  jet  of  2  C771 
in  radius  is  centered  at  x  =  2.75  cm.y  =  0  cm.  The  initial  temperature  of  the  jet  is  at 
1000  °K  at  1  atm. 
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Cliebyshev  collocation  methods  is  used  in  both  x  and  y  directions. 

The  Chebyshev  collocation  points  ^  is  mapped  into  another  set  of  interpolation  points 
X  by  a  grid  transformation  (Kosloff,  TahEzer)  in  the  form  of 


sin  HqQ 

sin“^(a:) 


Foi'  Of  =  sech((lne)/iV)  and  =  cos(7rj /N).  where  e  is  the  machine  zero.  CFL  number 
larger  than  the  standard  Chebyshev  methods  would  allow  can  be  taken.  For  N  =  512, 
CFL  can  be  15  times  larger  using  the  x{^j)  as  the  interpolation  points  without  any 
degradation  in  approximation  accuracy. 

Differentiation  of  the  three  dimensional  data  set  {n,m,k)  is  done  by  applying  Cosine 
Transform-Recursion-Inverse  Cosine  Transform  on  the  Cray  C90.  n,m  are  the  number 
of  points  in  x  and  y,  respectively,  k  is  the  number  of  PDE.  Typically,  n  —  m  =  512 
and  k  =  7. 

Boundary  Condition  : 

In  the  y  direction,  a  reflective  boundary  condition  is  used  at  y  =  ±6.5  cm.  Supersonic 
inflow  is  imposed  on  the  inflow.  At  the  outflow,  characteristic  treatment  based  on  the 
eigenfunctions  and  eigenvalues  of  the  one  dimensional  Euler  equation  is  used. 

A  low  storage  third  order  TVD  stable  Runge-Kutta  scheme  (Shu,  Osher)  is  used  to 
march  the  equations  in  time.  It  has  the  form  of 

=  W^-AtLiW^) 

+  -  AtL{W^) 

^  +  2W^  -2AtL{W^) 

where  L  is  the  spatial  operator  for  the  fluxes. 

At  each  Range  Kutta  time  step,  the  derivative  and  the  solution  are  filtered  by  the 
exponential  filter, 


a{k)  =  k  =  0,l,...,N 

where  a  =  —  In  e  and  e  is  the  machine  zero.  Typically,  7  =  14  for  the  smoothing  of  the 
derivative  of  the  flu.xes  and  7  =  12  for  the  smoothing  of  the  solution. 

It  is  mainly  used  for  the  purpose  of  stabilizing  the  spectral  scheme  from  Gibbs  phe¬ 
nomenon  and  nonlinear  instability. 

Gibbs  free  reconstruction  of  the  data  has  not  yet  been  implemented  in  a  satisfactory 
way.  Further  research  is  needed. 


•  At  some  stage  of  evolution  of  the  PDE,  1)  Gibbs  phenomena  and/or  small  scale  oscilla¬ 
tion  in  the  jets  caused  negative  hydrogen  mass  fraction,  2)  compression  of  the  jets  and 
lack  of  strong  numerical  dissipation  cause  unrealistic  local  large  gradients  in  species 
mass  fraction.  The  Global  nature  of  the  strong  smoothing  has  the  undesirable  effect 
of  destroying  information  in  regions  other  than  those  few  local  unstable  ones.  To  deal 
with  kind  of  situation,  points  that  exceed  a  certain  tolerance  level  are  smoothed  locally 
by  nine  point  averaging.  The  tolerance  is  set  by  limiting  the  hydrogen  mass  fraction 
to  not  exceed  one  percent  of  one  below  zero.  One  important  point  to  remember  is 
that  the  local  smoothing  is  done  on  the  primitive  valuables,  not  on  the  conservative 
valuables  (does  not  seem  to  work  well). 

1  Implementation  of  the  spectral  code 

•  The  Spectral  code  is  implemented  on  the  Gray  C90  at  WES  High  Performance  Gom- 
puting  (HPG)  Genter. 

•  The  code  has  about  88%  of  DO  loop  vectorized. 

•  The  code  make  used  of  the  multiple  fast  Fourier  Transform  routine  in  the  Cray  Scientific 
library  to  perform  the  main  computational  kernel,  namely,  the  differentiation  of  the 
fluxes  and  smoothing  of  the  solution. 

•  The  computational  kernel  of  the  algorithm  achieves  performance  of  550  megaflop  on  a 
single  processor  with  512  x  512  grids. 

•  Total  CPU  time  usage  for  the  512  x  512  grids  and  final  physical  time  T  =  80  microsec¬ 
ond  is  about  4|  CPU  hours.  Average  CPU  for  each  time  step  is  about  3.8  seconds. 
Total  number  of  time  step  is  about  4127. 

ENO  code 

•  Geometry  same  as  spectral,  except  y  G  [0  cm,  5  cm],  using  the  symmetry  property  of 
^he  flow. 

•  ENO  scheme  applied  to  a  regular  grid  with  mappings  applied  to  each  coordinate  di¬ 
rection  separately,  uniform  grid  on  x,  tM\{,8x)/c  in  y.  where  c  is  a  constant. 

•  Boundary  Conditions  : 

Reflective  boundaiy  conditions  at  y  =  0  cm  and  5  cm.  Characteristic  boundary  con¬ 
ditions  at  inflow  and  outflow. 

•  The  same  3rd  order  TVD  Runge-Kutta  applied  in  the  spectral  scheme  is  used  to 
advance  in  time,  using  the  method  of  lines. 

•  Each  chemical  species  is  treated  as  a  separate  conservative  variable  in  the  numerical 
formulation.  Negative  species  concentrations  are  treated  by  clipping  values. 
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•  The  hyperbolic  system  is  approximately  diagonalized  locally  in  order  to  decouple  the 
system  of  equations  in  the  ENO  derivative.  The  computed  Jacobian  matrix  includes 
contributions  due  to  the  chemical  species,  but  the  dependence  of  pressure  on  species 
mass  fraction  is  ignored. 

•  Local  Lax-Friedrich  flux-splitting  applied  for  upwinding. 

•  Viscosity  is  computed  by  central  differences.  If  s  is  the  shu  style  order  of  the  scheme, 
up  to  2'^(s+l)  order — user  chooses. 

3.1  ENO  Implementation 

•  The  ENO  code  is  implemented  on  the  Thinking  Machines  CM5  at  the  Army  High 
Performance  Computing  Research  Center. 

•  Code  implemented  using  the  CMMD  message  passing  library,  the  CDPEAC  assembly 
language  and  GCC-AC  vector  C  code. 

•  The  algorithm  achieves  in  excess  of  1  Gflop  on  a  32  node  machine  for  300  x  300  grids. 
On  a  512  node  machine,  the  code  achieves  in  excess  of  16  Gflops. 

•  CPU  time  usage  for  10000  time  steps  for  a  1472  x  736  grid  is  about  3  hours  for  a 
512  node  machine.  11029  time  steps  took  10243  seconds  (or  2.84  hours).  It  took  1.58 
gigabytes  of  memory,  total. 

4  Numerical  Results 

Here  we  would  like  to  study  the  resolution  properties  of  the  ENO  scheme  and  the  spectral 
scheme  for  supersonic  reactive  flows.  It  is  important  that  the  flame  front  (a  narrow  region 
where  the  combustion  take  place)  to  be  well  resolved.  The  physical  aspect  of  this  problem 
will  be  study  in  some  later  time.  Hence,  to  simplify  the  number  of  PDE  to  be  solved,  we 
solved  the  Euler  equation  without  combustion.  (ENO  version  of  the  full  combustion  model  is 
under  development.  The  spectral  version  of  the  full  combustion  model  is  ready  and  running.) 

The  2D  supersonic  reactive  flow  is  modeled  by  a  Mach  2  normal  air  shock  moving  toward 
right  side  of  the  domain  and  interacted  with  a  column  of  hjMrogen  jet  as  depicted  in  figure 
1.  (Solution  for  multiple  hydrogen  jets  configuration  is  also  available). 

The  hydrogen  mass  fraction  of  the  second  order  ENO  scheme  (Shu  s  definition)  with  grid 
size  240  x  120,  and  1472  X  736  are  shown  in  figures  2  and  3  respectively,  at  f  =  60  microsecond. 
The  corresponding  solution  from  the  spectral  code  with  grid  size  512  x  512  is  shown  in  figure 
5.  The  global  structure  (for  example,  density)  of  the  two  algorithms  is  very  similar.  The  fine 
structures  (for  example,  the  flame  front),  however,  was  not  captured  in  the  low  resolution 
version  of  the  ENO  scheme.  In  most  literature,  this  is  the  solution  that  other  researches 
(FCT,  ENO)  had  obtained.  As  the  resolution  of  the  ENO  scheme  increased,  it  then  able 
to  capture  the  finer  structure  inside  the  hydrogen  jet.  Still,  the  strong  inherent  numerical 
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dissipation  of  the  scheme  prevented  further  development  of  some  finer  structures  in  later 
time.  The  spectral  solution,  however,  shows  a  great  deal  more  structures.  More  noticeably 
is  the  mushroom  shape  structure  which  is  a  tell-tale  sign  of  Rayleigh-Taylor  instability. 

A  third  order  ENO  scheme  was  developed  with  mapping  to  cluster  points  near  the  jet. 
The  hydrogen  mass  fraction  now  indicate  a  completely  new  configuration  inside  the  hydrogen 
jet  than  the  low  resolution  case  (figure  5).  A  mushroom  shape  structure  is  now  emerged  from 
where  the  air  penetrating  the  hydrogen  (A  Rayleigh-Taylor  instability).  The  important  of 
high-order /high  resolution  scheme  in  resolving  the  combustion  front  is  evident  from  this 
study. 


2-fldadef  ENO  Scheme 
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